%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
% Use this file as a template for your Montesigliano report
% 
%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% 

\documentclass[graybox,envcountchap,sectrefs]{svmono}
\usepackage{mathptmx}
\usepackage{helvet}
\usepackage{courier}
\usepackage{type1cm}         
\usepackage{makeidx}         % allows index generation
\usepackage{graphicx}        % standard LaTeX graphics tool
\usepackage{multicol}        % used for the two-column index
\usepackage[bottom]{footmisc}% places footnotes at page bottom
\usepackage{import}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{algpseudocode} 
\usepackage{color}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\pgfplotsset{plot coordinates/math parser=false}
\usepackage{setspace,subfigure}
\usepackage{tikz}
\usetikzlibrary{plotmarks}
\usetikzlibrary{lindenmayersystems}
\usetikzlibrary{decorations.pathmorphing}
\newcommand{\comm}[1]{{\color{red}#1}}

\usepackage{algorithm}
\renewcommand{\algorithmicrequire}{\emph{function}}
\renewcommand{\algorithmicend}{\emph{\textbf{end}}}
\renewcommand{\algorithmicif}{\emph{\textbf{if}}}
\renewcommand{\algorithmicelse}{\emph{\textbf{else}}}
\renewcommand{\algorithmicfor}{\emph{\textbf{for}}}
\renewcommand{\algorithmicdo}{\emph{\textbf{do}}}
\renewcommand{\algorithmicwhile}{\emph{\textbf{while}}}
\renewcommand{\algorithmicthen}{\emph{\textbf{then}}}

% Colors
\def\red{\textcolor{red}} 
\def\blue{\textcolor{blue}}
\def\green{\textcolor{green}} 
\def\yellow{\textcolor{yellow}}

% parenthesis etc.
\def\f{\frac} 
\def\lb{\left[}  
  \def\rb{\right]} 
\def\lcb{\left\{} 
  \def\rcb{\right\}}
\def\lp{\left(}
  \def\rp{\right)} 
\def\la{\left\langle} 
  \def\ra{\right\rangle}

% hats etc.
\def\wt{\widetilde} 
\def\wh{\widehat} 

% definition
\def\v{\mathbf}
\newcommand{\bsigma}{\mbox{\boldmath$\sigma$}}
\newcommand{\bA}{{\bf A}}
\newcommand{\bS}{{\bf S}}
\newcommand{\bB}{{\bf B}}
\newcommand{\bC}{{\bf C}}
\newcommand{\bU}{{\bf U}}
\newcommand{\bP}{{\bf P}}
\newcommand{\bv}{{\bf v}}
\newcommand{\bc}{{\bf c}}
\newcommand{\bx}{{\bf x}}
\newcommand{\bI}{{\bf I}}


\begin{document}
\chapter{Fast methods for constructing scattering matrices for
  complicated scatterers in the plane}

\abstract{We describe an approach for solving a \red{scattering}
  problem in \red{the} plane. We consider scatterers that consist of a finite
  collections of hard inclusions. The developed method is capable of
  handling complicated aggregations of scatterers, such as scatterers
  that are within close proximity of each other. Our approach is based
  on fast direct methods for constructing scattering matrices and
  multilevel representation of operators. Our method achieves high
  accuracy and is applicable to large-scale problems. We illustrate
  the efficiency of our method by examin\red{ing} several numerical
  experiments that include various shapes and numbers of scatterers in
  a plane.}

\section{Introduction}

This work has been a result of a project initiated in the
Montestigliano workshop on \emph{fast methods for scientific
  computing}~\cite{montestigliano12}. In the spring of 2012 the topic
of the workshop focuses on developing fast solvers for elliptic
partial differential equations such as the Laplace and Helmholtz
equations, which are known to form some of the key building blocks in
scientific computing. The workshop introduces a series of linear
complexity direct solvers that can under certain conditions outperform
state-of-the-art iterative methods. These methods are robust, stable,
and fast. In the present work, we utilize fast direct methods for
constructing scattering matrices to solve a scattering problem in a
plane with complicated domain formulation.

We consider here a class of constant coefficient Helmholtz equation
given on a planar \red{unbounded} domain. In addition, we concentrate
on scatterers with domains that are smooth, i.e. without any corners
(for an extension to the case with corners see~\cite{bre11}). For this
set of problems, it is advantageous to reformulate the boundary value
problem into its equivalent integral equation and our method is based
on such an approach. One of the major advantages of this approach is
that the integral operator given on a Lipschitz curve is
well-conditioned when the underlying function space is
squared-integrable~\cite{bre11}. However, the transformation of a
boundary value problem to its integral formulation does not
immediately provides an computationally efficient method for its
solution. The discretization of the corresponding integral equation
usually yields a dense system of linear equations and in majority of
the cases, the solution is computationally expensive. Traditional
methods utilize iterative solvers, e.g. GMRES and multigrid, to
improve the computational time~\cite{rokgre87,rok90,rok93}, but they
still fall short of what is required for large-scale problems.

Recent developments have enabled the construction of direct solvers
that can accurately obtain solution to large-scale system of linear
equations with linear complexity. These methods, called fast direct
solvers, can construct an approximation to the solution operator of
the problem in a single step. The present method utilizes the fast
direct methods described in~\cite{marrok05,gilyoumar11} to construct
\textit{scattering matrices} for solving the scattering problem.

The utility of scattering operators provides a mechanism for
accelerating the numerical calculations necessary to solve
non-oscillatory elliptic boundary value problems. The concept of a
scattering operator is very powerful and is used in acoustics and
electromagnetics to great effect. Here, we show that the scattering
matrices exhibit desirable properties for computations. In particular,
we demonstrate that the scatter\red{ing} matrices can often be
represented as structured matrices. Furthermore, the proposed method
use domain decomposition techniques to divide the computational domain
into a number of small patches. On each patch, a so called \emph{proxy
  matrix} is constructed. This proxy matrix contains all information
about the interior geometry of the patch that is needed to solve the
global system. The proxy matrix can be thought of as a compact
representation of the scattering matrix of the patch.

The rest of the paper is organized as follows: in
Section~\ref{sec.problem-formulation}, we formulate the scattering
problem and describe the numerical method used to solve it; in
Section~\ref{sec.results}, we demonstrate the utility of our method on
several numerical experiments; and finally in
Section~\ref{sec.conclusion}, we summarize our method and results, and
discuss possible future directions.

\section{Problem Formulation}
\label{sec.problem-formulation}

Consider the following boundary value problem
\begin{subequations}
  \label{eq.helmholtz-bc}
  \begin{align}
    \label{eq.helmholtz}	
    -\Delta \, u(\mathbf{x}) \, - \, k^{2} \, u(\mathbf{x})
    & 
    \;\; = \;\;
    f(\mathbf{x}),
    &
    \mathbf{x} \, \in \, \Omega,
    \\[0.05cm]
    \label{eq.bc}
    u(\mathbf{x})
    & 
    \;\; = \;\;
    0,
    &
    \mathbf{x} \, \in \, \Gamma,
  \end{align}
\end{subequations}
where $\Omega$ is an unbounded domain that is exterior to a smooth
contour $\Gamma$, $u$ is the total wave field, $f$ is the body load,
and $k$ is the wavenumber; see figure~\ref{fig.scattering_problem} for
an illustration. Let $v$ be a field generated by $f$ which satisfies
the following boundary value problem,
\begin{equation}
  -\Delta v(\mathbf{x}) \, - \, k^2 \, v(\mathbf{x}) \; = \; f(\mathbf{x}), \;\;\; \mathbf{x} \in \mathbb{R}^2.
  \label{eq.incoming-field}
\end{equation}
Physically, $v$ corresponds to a wave incoming on $\Gamma$. A
scattered wave needs to be generated to satisfy the boundary
condition~\eqref{eq.bc}. Let the total field, $u$, be a summation of the
incoming, $v$, and outgoing, $w$, fields
\begin{equation}
  u(\mathbf{x}) \; = \; v(\mathbf{x}) \, + \, w(\mathbf{x}),
  \label{eq.total-field}
\end{equation}
then the radiated wave $w$ solves the following boundary value problem,
\begin{subequations}
  \label{eq.outgoing-field-bc}
  \begin{align}
    \label{eq.outgoing-field}	
    -\Delta \, w(\mathbf{x}) \, - \, k^{2} \, w(\mathbf{x})
    & 
    \;\; = \;\;
    0,
    &
    \mathbf{x} \, \in \, \Omega,
    \\[0.05cm]
    \label{eq.w-bc}
    w(\mathbf{x})
    & 
    \;\; = \;\;
    -v(\mathbf{x}),
    &
    \mathbf{x} \, \in \, \Gamma.
  \end{align}
\end{subequations}
For well-posedness, we imposed a radiation condition at infinity on \eqref{eq.outgoing-field-bc},
\begin{equation}\label{eq.rad_cond}
  \lim_{|\bx| \to \infty} \sqrt{|\bx|} \left(
    \frac{\partial}{\partial |\bx|} \, - \, \mathrm{i} k \right) w(\bx) \; = \; 0,
\end{equation}
which indicates that the field is outgoing.

\begin{figure}[t]
  \begin{center}
    \begin{tabular}{cc}
      \subfigure[]
      {
        \includegraphics[width=0.4\textwidth,keepaspectratio]{figures/scattering/figscattering1.pdf}
        \label{fig.scattering_problem}
      }
      &
      \subfigure[Incoming field $v$]
      {
        \includegraphics[width=0.4\textwidth,keepaspectratio]{figures/inouttot/figinouttotin.pdf}
        \label{fig.scatteringin}
      }\\
      \subfigure[Outgoing field $w$]
      {
        \includegraphics[width=0.4\textwidth,keepaspectratio]{figures/inouttot/figinouttotout.pdf}
        \label{fig.scatteringout}
      }
      &
      \subfigure[Total field $u=v+w$]
      {
        \includegraphics[width=0.4\textwidth,keepaspectratio]{figures/inouttot/figinouttottot.pdf}
        \label{fig.scatteringtot}
      }
    \end{tabular}
  \end{center}
  \vspace{-0.3cm}
  \caption{(a) A smooth contour $\Gamma$ is enclosed in an unbounded
    domain $\Omega$. (b) An incoming field, $v$, generated by external
    sources approaches $\Gamma$ which induces charges on $\Gamma$. (c)
    The induced charges on $\Gamma$ generate an outgoing field
    $w$. (d) The total field is a summation of the incoming and
    outgoing fields. Figures are generated by the {\sc Matlab} file
    \texttt{ex\_inouttot(3,5*pi/4)}.}
  \label{fig.scattering}
\end{figure}

\subsection{Numerical Formulation}
\label{sec.numerical-formulation}

\begin{figure}[t]
  \begin{center}
    \includegraphics[width=0.35\textwidth,keepaspectratio]{figures/scattering/figscattering2.pdf}
  \end{center}
  \vspace{-0.3cm}
  \caption{A smooth contour $\Gamma$ is enclosed in an unbounded
    domain $\Omega$. The matrix $\bC$ maps a charge distribution at
    the points on $D$ to a potential at the discretization points on
    $\Gamma$. The matrix $\bB$ maps the charges at points on $\Gamma$
    to a potential at points on $D$. The scattering matrix $\bS$ is
    loosely defined as the matrix that maps the incoming field $v$ to
    the outgoing field $w$.}
  \label{fig.single_contour}
\end{figure}

We represent the outgoing field $w$ using a layer potential
\begin{equation}
  w(\mathbf{x}) \; = \; \int_{\Gamma} K(\mathbf{x},\mathbf{y}) \, \sigma(\mathbf{y}) \mathrm{d}s(\mathbf{y}), \;\;\; \mathbf{x} \, \in \, \Omega,
\end{equation}
where $K$ denotes the kernel given by
\begin{equation}
  K(\mathbf{x},\mathbf{y}) \; = \; {\bf n}(\mathbf{y}) \cdot \nabla_{\mathbf{y}} \, \phi(\mathbf{x} - \mathbf{y}) \, + \, \mathrm{i} k \, \phi(\mathbf{x} - \mathbf{y}),
\end{equation}
the vector ${\bf n}(\mathbf{y})$ is the outwards pointing unit
normal to $\Gamma$ for $\mathbf{y} \in \Gamma$. The source
distribution $\sigma$ satisfies the second kind Fredholm equation
\begin{equation}
  2 \pi \mathrm{i} \sigma(\mathbf{x}) \, + \, \int_{\Gamma} K(\mathbf{x},\mathbf{y}) \, \sigma(\mathbf{y}) \mathrm{d}s(\mathbf{y}) \; = \; -v(\mathbf{x}), \;\;\; \mathbf{x} \, \in \, \Gamma.
  \label{eq.BIE}
\end{equation}

In order to discretize the boundary integral equation, let $\{ \mathbf{x}_{i},
m_{i} \}_{i = 1}^{N}$ denote a quadrature rule on $\Gamma$ so that for
any smooth $\varphi$ we have
\begin{equation}
  \int_{\Gamma} \varphi(\mathbf{x}) \mathrm{d}s(\mathbf{x}) \; \approx \; \sum_{i = 1}^{N} m_{i} \varphi(\mathbf{x}_{i}).
\end{equation}
Then using Nystr$\ddot{\text{o}}$m discretization, we obtain the
discretized equation
\begin{equation}
  {\bA \bsigma} \; = \; - {\bv},
\end{equation}
where ${\bf v}$ is the vector given by
$
\bv_i \; = \; v(\mathbf{x}_{i}),
$
and $\bsigma$ is a vector that approximates the exact solution
$
\bsigma_i \; \approx \; \sigma(\mathbf{x}_{i}).
$
Here ${\bA}$ is an $N \times N$ matrix with entries
\begin{equation}
  \bA_{i j} \; = \; 2 \pi \mathrm{i} \delta_{i j} \, + \, m_j \, K(\mathbf{x}_{i}, \mathbf{x}_{j}),
\end{equation}
when $\mathbf{x}_{i}$ and $\mathbf{x}_j$ are not close. When
$\mathbf{x}_{i}$ and $\mathbf{x}_j$ are close, the entries of $\bA$
must be modified because the kernel $K$ is singular.


\subsection{Scattering Operators}
\label{sec.scattering-matrices}

Scattering operators have been successfully used as a tool for
accelerating numerical calculations~\cite{bre11}. The scattering
operator $S$ associated with the contour $\Gamma$ is loosely defined
as the operator that maps the incoming field, $v$, to the outgoing
field, $w$. To formalize the definition, let $D$ denote a disc that is
large enough to hold the scatterer $\Gamma$ and is small enough that
$f$ is supported outside of $D$; see figure~\ref{fig.single_contour}
for an illustration.

Let $\phi$ represent the fundamental solution of the
Helmholtz operator~(\ref{eq.helmholtz}):
\begin{equation}
  \phi(\mathbf{x}) \; = \; \cfrac{\mathrm{i}}{4} \; H_{0}^{(1)} \left( k |\mathbf{x}| \right),
  \label{eq.sol-helmholtz}
\end{equation}
where $\mathrm{i}$ is the imaginary unit and $H_{0}^{(1)}$ is the
Hankel function of the first kind of order $0$. The fundamental
solution~(\ref{eq.sol-helmholtz}) satisfies the radiation condition at
infinity \eqref{eq.rad_cond}, and the partial differential equation
\begin{equation}
  -\Delta \phi(\mathbf{x}) \, - \, k^2 \, \phi(\mathbf{x}) \; = \; \delta(\mathbf{x}),
\end{equation}
where $\delta(\mathbf{x})$ is the Dirac function. The incoming field
can be represented using a charge distribution $q$ on the boundary of
the disk $\partial D$,
\begin{equation}
  {\displaystyle v(\mathbf{x}) \; = \; \int_{\partial D} \phi (\mathbf{x} - \mathbf{y}) \, q(\mathbf{y}) \, \mathrm{d}s(\mathbf{y})}.
\end{equation}

% The incoming field defined by~(\ref{eq.incoming-field}) can be
% determined by
% \begin{equation*}
%   v(\mathbf{x}) \; = \; \left[ v \circ f \right](\mathbf{x}) \; = \; \int_{\mathbb{R}^{2}} \phi (\mathbf{x} - \mathbf{y}) \, f(\mathbf{y}) \, \mathrm{d}\mathbf{y},
% \end{equation*}
% where $\circ$ denotes the convolution operator. In addition, since $f$
% is supported outside of the domain $D$,
% equation~(\ref{eq.outgoing-field}) converges to
% \begin{equation*}
%   -\Delta v(\mathbf{x}) \, - \, k^2 \, v(\mathbf{x}) \; = \; 0, \;\;\; \mathbf{x} \, \in \, D.
% \end{equation*}
% It is worthy to note that the incoming field $v$ represents the
% particular solution to~(\ref{eq.helmholtz}) that accounts for the body
% load $f$. However, it does not satisfy the homogeneous Dirichlet
% boundary conditions on $\Gamma$. The outgoing field $w$ can be thought
% of as a correction that ensures that the total field $u$ satisfies the
% boundary condition.

% The outgoing field $w$ is represented simply by tabulation on
% $\partial D$. 
The mapping from the incoming to the outgoing field can be represented
in terms of the scattering operator $S$ defined as
\begin{equation}
  S \, : \, q|_{\partial D} \; \mapsto \; w|_{\partial D}.
\end{equation}

We note that the representations of the incoming and outgoing fields
using charge distribution on $\partial D$ and tabulation on $\partial
D$ are complicated by the fact that the Helmholtz problem can be
resonant on $D$. For brevity, we will not address this issue here.

Let $\{ \v z_{i}, m_{i} \}_{i = 1}^{M}$ denote a set of
discretization points and a quadrature rule on $\partial D$. In a
discrete form, the scattering matrix can be written as
\begin{equation}
  \bS \; = \; \bB \, \bA^{-1} \, \bC,
  \label{eq.Smatrix}
\end{equation}
where $\bA$ is the discretized boundary integral
operator~(\ref{eq.BIE}), $\bC$ is the matrix that maps a charge
distribution $\v q$ at points $\{\v z_i \}_{i = 1}^{M}$ to a potential
at the discretization points $\{\v x_i \}_{i = 1}^{N}$, and $\bB$ is
the matrix that maps the charges $\v \sigma$ at points $\{\v x_i \}_{i
  = 1}^{N}$ to a potential at points $\{\v z_i \}_{i = 1}^{M}$. The
matrices $\bC$ and $\bB$ are given by
\begin{equation}
  \bC_{ij} \; = \; \cfrac{2 \pi}{M} \, \phi(\mathbf{x}_{i}, \mathbf{z}_{j}),
  \;\;\;\;\;
  \bB_{ij} \; = \; m_j \, K(\mathbf{z}_i, \mathbf{x}_j),
\end{equation}

\subsection{Proxy Matrices}
\label{sec.proxy-matrices}

We note that the matrices $\bC$ and $\bB$ are typically
rank-deficient. This means that they each admits an interpolatory
decomposition~\cite{marrok05}
\begin{equation}
  \begin{array}{ccccccccc}
    \bC
    & = &
    \bU
    &
    \tilde{\bC},
    & \;\;\;\;\;\; &
    \bB
    & = &
    \tilde{\bB}
    &
    \bU^{*},
    \\[0.2cm]
    \scriptscriptstyle N \times M
    & &
    \scriptscriptstyle N \times s
    &
    \scriptscriptstyle s \times M
    & \;\;\;\;\;\; &
    \scriptscriptstyle M \times N
    & &
    \scriptscriptstyle M \times s
    &
    \scriptscriptstyle s \times N
  \end{array}
  \label{eq.CB-id}
\end{equation}
where $s$ is the numerical rank of matrices $\bC$ and $\bB$, $\bU$ is
the base matrix that contains a $s \times s$ identity matrix and no
entry is larger than one, and matrices $\tilde{\bC}$ and $\tilde{\bB}$ are the
skeletons of $\bC$ and $\bB$, i.e., $\tilde{\bC} = \bC(\bI_s, :)$ and
$\tilde{\bB} = \bB(:, \bI_s)$ with $\bI_s$ denoting an index vector that
identifies the skeleton subset of the discretization points $\{ \bx_{i}
\}_{i=1}^{N}$ on $\Gamma$.

Substitution of~(\ref{eq.CB-id}) for $\bC$ and $\bB$
into~(\ref{eq.Smatrix}) yields
\begin{equation}
  \bS \; = \; \tilde{\bB} \, \bU^{*} \, \bA^{-1} \, \bU \, \tilde{\bC} \; = \;
  \tilde{\bB} \, \bP \, \tilde{\bC},
\end{equation}
where $\bP$ is defined as the proxy matrix for $\Gamma$. For practical
purposes, the proxy matrix $\bP$ contains all the information that is
needed from the scattering matrix $\bS$. The flanking matrices
$\tilde{\bC}$ and $\tilde{\bB}$ are essentially just clutter.


\section{Results}
\label{sec.results}

In this section, we illustrate the utility of our method on several
numerical experiments that include various shapes and numbers of
scatterers. We first discuss the \red{various?} parameters (e.g.,
\red{wavenumber}, number of discrete points, \red{shape} of scatterer)
that can affect the numerical rank of the scattering matrix. Then, a
divide-and-conquer method is introduced to solve problems with
multiple scatterers within close proximity of each other.

% \subsection{\comm{Low Rank Approximation of the Scattering Matrix Valid in
% the Far Field, skip since it comes in later numerically?}}

\subsection{Parameters Affecting the Numerical Rank of the Scattering Matrix}

Here, we analyze how various parameters can affect the number of nodes
that are retained in the skeleton of the original scattering
matrix. Parameters of interest include the size of the contour $D$
(\red{in this section, a} disc with radius $R$) and the wavenumber
$k$. In addition, we study various shapes of the scatterer that can
impact the numerical rank.

We first \red{set} the wavenumber \red{to} $k = 1$ and study the
effect of contour\red{s} $D$ of various sizes and four different
scatterers on the the numerical rank of the scattering
matrix. Figure~\ref{fig:skeleton}(a) shows the number of nodes that
are retained in the skeleton of the original scattering matrix for
various sizes of the disc $D$ and different shapes of the
scatterer. At approximately $R = 3R_{\mathrm{min}}$ we attain the
minimal number of nodes required in the skeleton for a given accuracy
$\varepsilon$, as depicted in figure~\ref{fig:skeleton}(a). Here,
$R_{\mathrm{min}}$ denotes the minimum radius of a circle $D$ that
contain\red{s} the scatterer.  Increasing the radius of the contour
$D$ \red{further} does not decrease the number of skeleton nodes.
This optimal radius (since a minimal number of nodes are required on
the disc $D$) decreases as $k$ increases, as shown in
figure~\ref{fig:skeleton}(b).

Furthermore, as the contour of the scatterer becomes more complex,
more skeleton nodes are required to attain a specified numerical
tolerance. However, the difference in the number of skeleton nodes are
minor between the different shapes.

% The formulas for the investigated scatterers are given below, star, moon and an ellipse:
% \[
% \mathbf{x}_{\mathrm{star}}=\mathbf{x}_c+
% \left(\begin{array}{l}
%     1.5\cos(t)+\frac{r}{2}\cos((p+1)t)+\frac{r}{2}\cos((p-1)t)\\
%     \sin(t)+\frac{r}{2}\sin((p+1)t)-\frac{r}{2}\sin((p-1)t)
%   \end{array}\right)
% \]
% \[
% \mathbf{x}_{\mathrm{moon}}=\mathbf{x}_c+
% \left(\begin{array}{l}
%     0.5\cos(t)+\cos^2(t)\\
%     \sin(t)
%   \end{array}\right)
% \quad
% \mathbf{x}_{\mathrm{ellipse}}=\mathbf{x}_c+
% \left(\begin{array}{l}
%     \alpha\cos(t)\\
%     \beta\sin(t)
%   \end{array}\right)
% \]
% \[
% t=0,\ldots,2\pi,
% \quad
% \mathbf{x}_c=(0,0)^{\mathrm{T}},\quad
% r=0.3,\quad p=\{5,10\},
% \]
% where $\bx$ is the center of the object, \comm{$r$ controls the depth of the ``arms'' of the star}, $p$ is the number of points of the star, $\alpha$ and $\beta$ are the semi-major and semi-minor axis of the ellipse.
\begin{figure}[t]
  \centering
  \begin{tabular}{cc}
    \subfigure[]{\includegraphics[width=0.45\textwidth,keepaspectratio]{figures/skeleton/figskelr.pdf}}& \subfigure[]{\includegraphics[width=0.45\textwidth,keepaspectratio]{figures/skeleton/figskelk.pdf}}\\
  \end{tabular}
  \vspace{-0.2cm}
  \caption{(a): Four different objects, $(M=500, N=400, k=1,
    \varepsilon=10^{-10})$, Varying radius $R$ of \red{the} disc $D$,
    depending on $R_{\mathrm{min}}$. (b): Five pointed star, $(M=500,
    N=400$, $k=1,5,10,\ldots,45, \varepsilon=10^{-10})$. As $k$
    increases the number of nodes of the skeleton increases. Varying
    radius $R$ of disc $D$, depending on $R_{\mathrm{min}}$. Figures
    generated with \texttt{ex\_skeleton}.}
  \label{fig:skeleton}
\end{figure}

We next study the effects of different wavenumbers $k$ for a single
five-pointed smooth star scatterer and disc
$D$. Figure~\ref{fig:skeleton}(b) shows the number of skeleton nodes
retained in the original scattering matrix for several values of
$k$. Note that as $k$ increases the number of skeleton nodes increases
for a fixed $\varepsilon$. The increase in the \red{number of}
skeleton nodes exhibits a linear dependence to the increase in $k$. In
addition, the optimal radius $R$ of $D$ behaves linear in $k$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Two Scatterers: Merging Skeletons}

We have shown in the previous section that the mapping between a
source distribution outside of $D$ to the wave scattered by $\Gamma$
can be accurately represented using only the skeleton points. Here, we
consider the situation when we have two scatterers, $\Gamma_1$ and
$\Gamma_2$, as depicted in figure~\ref{fig:2obj}.
\begin{figure}[t]
  \centering %\import{./figures/}{merge.pdf_tex}
  \includegraphics[width=0.5\textwidth,keepaspectratio]{figures/2obj/fig2obj.pdf}
  \vspace{-0.2cm}
  \caption{Scattering by two objects.}
  \label{fig:2obj}
\end{figure}

A proxy matrix for the matrix $\bA$ can be obtained by treating both
scatterers as a single scatterer, $\Gamma_{\mathrm{tot}} = \Gamma_1
\cup \Gamma_2$. Let $N$ denote the number of discretization points of
each scatterer, then the computation of the proxy matrix for $\bA$
would require computing an inverse and an interpolatory decomposition
of a matrix of size $2N \times 2N$. However, the computational cost
would be approximately eight times that associated with a single
object (the inversion of a dense $M\times M$ matrix requires order
$N^3$ operations). Instead, we propose to utilize the proxies computed
separately for $\Gamma_1$ and $\Gamma_2$ to build a proxy for
$\Gamma_{\mathrm{tot}}$. The discrete problem for two scatterers can
be decomposed into a block form
\begin{equation} 
  \label{eq:2obj}
  \left[
    \begin{array}{cc} \bA_{11} & \bA_{12} \\ \bA_{21} & \bA_{22}\end{array}
  \right] 
  \left[
    \begin{array}{c} \sigma_1 \\ \sigma_2 \end{array}
  \right]  
  \; = \;
  - \left[ 
    \begin{array}{c} \bv_1 \\ \bv_2 \end{array}
  \right],
\end{equation}
where the blocks $\bA_{ii}$ correspond to the discretization of the
scattering problems with only $\Gamma_i$, and $\bA_{ij}$ map a source
distribution on $\Gamma_{j}$ to a potential on $\Gamma_i$. Let $D_1$
be a contour containing $\Gamma_1$ such that $\Gamma_2$ is outside of
$D_1$ and $D_2$ be a contour containing $\Gamma_2$ such that
$\Gamma_1$ is outside of $D_2$; see figure~\ref{fig:2obj} for an
illustration. Then utilizing the skeleton points to capture the effect
of $\Gamma_1$ and $\Gamma_2$ on the source distributions $\sigma_1$
and $\sigma_2$, respectively,~\eqref{eq:2obj} can be represented by
\begin{equation} 
  \label{eq:2obj-proxies}
  \underbrace{
    \left[ 
      \begin{array}{cc} \bP_{1}^{-1} & \bA_{12}(\bI_1,\bI_2) 
  	\\  \bA_{21}(\bI_2,\bI_1) & \bP_{2}^{-1}\end{array}
    \right] 
  }_{\tilde{A}}
  \underbrace{
    \left[
      \begin{array}{c} \sigma_1(\bI_1) \\ \sigma_2(\bI_2) \end{array}
    \right]
  }_{\tilde{\sigma}}  
  \; = \;
  -\underbrace{
    \left[ \begin{array}{c} \bv_1(\bI_1) \\ \bv_2(\bI_2) \end{array}
    \right]}_{\tilde{v}},
\end{equation}
where $\bP_i$ is the proxy matrix for $\Gamma_i$ and $D_i$, and
$\bI_i$ is the vector containing the indices of the corresponding
skeleton points.

Figure~\ref{fig:2obj_err}$(a)$ shows that \red{when the above method
  is used} the error \red{on the scattered field} is smaller than the
required accuracy, $\varepsilon=10^{-10}$, outside of the two
contours. Using the union of the skeletons for $\Gamma_1$ and
$\Gamma_2$, a new skeleton can then be built for
$\Gamma_{\mathrm{tot}}$ (based on an outer contour $D_{\mathrm{tot}}$)
such that the charge distribution at these points
$\tilde{\tilde{\sigma}}$ satisfies
\begin{equation}\label{eq:2obj-proxy}
  \tilde{\tilde{\bA}} \tilde{\tilde{\sigma}} = -\tilde{\tilde{\bv}}.
\end{equation}
\red{If $D_{\mathrm{tot}}$ contains $D_1$ and $D_2$ then the latter
  approximation allows an accurate computation of the radiated field
  outside of $D_{\mathrm{tot}}$, as shown in
  figure~\ref{fig:2obj_err}$(b)$.}


% Equation \eqref{eq:2obj} can be rewritten as
% \begin{align} 
%   A_{11} \sigma_1 & = -v_1 - A_{12} \sigma_2 \label{eq:2obj-1}
%   \\
%   A_{22} \sigma_2 & = -v_2 - A_{21} \sigma_1 \label{eq:2obj-2}
% \end{align}
% Let $D_1$ be a contour containing $\Gamma_1$ such that $\Gamma_2$ is
% outside of $D_1$ and $D_2$ be a contour containing $\Gamma_2$ such
% that $\Gamma_1$ is outside of $D_2$; see figure~\ref{fig:2obj} for an illustration. Let $P_i$ be the proxy matrix for $\Gamma_i$ and $D_i$, and let $I_i$ be the
% indices of the corresponding skeleton points.
% 
% The sources on the RHS of \eqref{eq:2obj-1} are located outside of
% $D_1$, and the potential scattered by $\Gamma_1$ in \eqref{eq:2obj-1}
% (the term $A_{21} \sigma_1$) is considered outside of
% $D_1$. Consequently in order to capture the effect of $\Gamma_1$ on
% the source distribution $\sigma_2$ on $\Gamma_2$ can be accurately
% represented using only the skeleton points on $\Gamma_1$. The same
% thing applies to the effect of $\Gamma_2$ on the scattering by
% $\Gamma_{\mathrm{tot}}$, so \eqref{eq:2obj} can be replaced by 
% \begin{equation} \label{eq:2obj-proxies}
%   \underbrace{\lp \begin{array}{c|c} P_{1}^{-1} & A_{12}(I_1,I_2) \\ \hline A_{21}(I_2,I_1) & P_{2}^{-1}\end{array}
%   \rp }_{\tilde{A}}
%   \underbrace{\lp \begin{array}{c} \sigma_1(I_1) \\ \sigma_2(I_2) \end{array}
%   \rp}_{\tilde{\sigma}}  = -
%   \underbrace{\lp \begin{array}{c} v_1(I_1) \\ v_2(I_2) \end{array}
%   \rp}_{\tilde{v}}
% \end{equation}
% in order to evaluate the scattered field outside of $D_1\cup
% D_2$. 

\begin{figure}[]
  \centering
  \begin{tabular}{cc}
    \subfigure[$\log_{10}|\phi_\sigma - \phi_{\tilde{\sigma}}|$]{\includegraphics[width=.45\textwidth,keepaspectratio]{figures/2objerr/fig2objerr.pdf}}&
    \subfigure[$\log_{10}|\phi_\sigma - \phi_{\tilde{\tilde{\sigma}}}|$]{\includegraphics[width=.45\textwidth,keepaspectratio]{figures/2objerr/fig2objerr2.pdf}}\\
  \end{tabular}
  \vspace{-0.15cm}
  \caption{Approximation errors. $\phi_{\sigma}$ is computed using the
    charge distribution $\sigma$ on all the discretization points of
    $\Gamma_{\mathrm{tot}}$ given by \eqref{eq:2obj}.  $\phi_{\tilde{\sigma}}$ is
    computed using the charge distribution $\tilde{\sigma}$ on the skeleton
    points of $\Gamma_1$ and $\Gamma_2$, given by
    \eqref{eq:2obj-proxies}. $\phi_{\tilde{\tilde{\sigma}}}$ is computed using
    the charge distribution $\tilde{\tilde{\sigma}}$ on the skeleton points of
    $\Gamma_{\mathrm{tot}}$, given by \eqref{eq:2obj-proxy}. $(a)$:
    $\log_{10}|\phi_\sigma - \phi_{\tilde{\sigma}}|$. The red rectangles
    correspond to the contours $D_1$ and $D_2$. $(b)$:
    $\log_{10}|\phi_\sigma - \phi_{\tilde{\tilde{\sigma}}}|$. The red circle
    correspond to the contour $D_{\mathrm{tot}}$. Figures are generated by
    the {\sc Matlab} file \texttt{twoobj(1.,0.)}.}
  \label{fig:2obj_err}
\end{figure}

\subsection{Multiple Scatterers: Hierarchical Decomposition}

In this section, we extend the procedure described in the previous
section for creating the scattering matrix for two objects to an
arbitrary number of scatterers. For simplicity, let us assume that
there are \red{$N_{obj} = 2^n$ objects. If a skeleton is build for
  each object, and these skeleton are then grouped two by two, then
  one is left with $2^{n-1}$ objects. The same idea can then be
  applied again and again until only one set of skeleton points are
  left. The algorithm is given
  algorithm~\ref{alg.multiple-scatterers}. In the pseudo-code $C_i$
  denotes the initial discretization points on the $i^{th}$ contour
  $\Gamma_i$.}

\begin{algorithm}
  \caption{Pseudo-code for building the proxy matrix for $2^n$
    scatterers.}
  \label{alg.multiple-scatterers}
  \begin{algorithmic}[1]
    \For{$i = 1,...,N_{\mathrm{obj}}$}
    \State Build the $\bA_{ii}$ 
    \EndFor
    
    \While{$N_{\mathrm{obj}} > 1 $}
    \For{$i = 1,3,...,N_{\mathrm{obj}}-1$}
    \State Build a contour $D_i$
    \State Build the proxy $\bP_{i}$ (Skeleton $\tilde{C}_i$)
    \State Build a contour $D_{i+1}$
    \State Build the proxy $\bP_{i+1}$ (Skeleton $\tilde{C}_{i+1}$)
    \State Build a contour $D_{i+1}$
    \State Get $\tilde{\bA}_{i,i+1}$
    \EndFor
    \State $N_{\mathrm{obj}} = N_{\mathrm{obj}}/2$
    \For{$i = 1,..,N_{\mathrm{obj}}/2$}
    \State $\bA_i = \tilde{\bA}_{2i,2i+1}$
    \State $C_i =\tilde{C}_{2i} \cup \tilde{C}_{2i+1}$
    \EndFor
    \EndWhile
    \State Build a contour $D_0$
    \State Build the proxy $\bP_0$ (Skeleton $\tilde{C}_0$)
  \end{algorithmic}
\end{algorithm}

For this procedure to produce a valid approximation, it is necessary
that every time a contour $D_i$ is chosen all of the points in $C_{j
  \ne i}$ are outside of it. Figure~\ref{fig:4obj} presents a
situation where four objects are located such that a rectangle can be
chosen for $D$ at each iteration.

% \begin{figure}[]
%   \begin{center}
%     \begin{tabular}{cc}
%       \includegraphics[width=.45\textwidth,keepaspectratio]{./figures/8obj_level01}&
%       \includegraphics[width=.45\textwidth,keepaspectratio]{./figures/8obj_level02}\\
%       $(a)$ : initial discretization & $(b)$: level 1\\
%       \includegraphics[width=.45\textwidth,keepaspectratio]{./figures/8obj_level03}&
%       \includegraphics[width=.45\textwidth,keepaspectratio]{./figures/8obj_level04}\\
%       $(c)$ : level 2 & $(d)$: level 3\\
%     \end{tabular}
%     \includegraphics[width=.45\textwidth,keepaspectratio]{./figures/8obj_level05}\\
%     $(e)$ : level 4
%   \end{center}
%   
%   \caption{Evolution of the number of discretization points at each
%   iteration of the hierarchical procedure. Figures generated by
%   \texttt{eightobj(1.,0.)}.\comm{do this for 4 objects to get 4 figs instead of 5?}}
%   \label{fig:8obj}
% \end{figure}

%%%%%%% alternative version of 4 scatterer visualization %%%%%%%%
\begin{figure}
  \begin{center}
    \begin{tabular}{cc}
      \subfigure[Initial discretization]{\includegraphics[width=.4\textwidth,keepaspectratio]{figures/4obj/fig4obj1.pdf}}&
      \subfigure[Level 1]{\includegraphics[width=.4\textwidth,keepaspectratio]{figures/4obj/fig4obj2.pdf}}\\
      \subfigure[Level 2]{\includegraphics[width=.4\textwidth,keepaspectratio]{figures/4obj/fig4obj3.pdf}}&
      \subfigure[Level 3]{\includegraphics[width=.4\textwidth,keepaspectratio]{figures/4obj/fig4obj4.pdf}}\\
    \end{tabular}
  \end{center}
  \caption{Evolution of the number of discretization points at each
    iteration of the hierarchical procedure. Figures are generated by the {\sc Matlab} file
    \texttt{ex\_fourobj(1.,0.)}.}
  \label{fig:4obj}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Closely Situated Scatterers}

Every time a proxy is build for an object $C_i$, an external contour
$D_i$ needs to be set such that all of the points in $C_{j \ne i}$ are
outside of $D_i$. In the example of figure~\ref{fig:4obj} objects are
chosen such that a rectangle can be chosen for each of the
$D_i$'s. This is in general not possible 

\red{In general it is not possible to find rectangles $D_i$ containing
  the contours $\Gamma_i$ that are such that all of the points in
  $C_{j \ne i}$ are outside of $D_i$ (see
  figure~\ref{fig:extcontour})}. Any shape is possible for $D_i$, but
a complex shape such as in figure~\ref{fig:extcontour}$(a)$ would not
be convenient. Instead, one can first choose a rectangle $\tilde{D}$,
and then choose $D_i$ to be the union of $\tilde{D}$ and of all the
points of $C_{j \ne i}$ interior to $D_i$ (situation represented in
figure~\ref{fig:extcontour}$(b)$). This way some of the points on
$D_i$ will not be useful (top and bottom right corners in
figure~\ref{fig:extcontour}$(b)$) but the effect of object $i$ onto
the other objects will be correctly taken into account.

\begin{figure}
  \begin{center}
    \begin{tabular}{cc}
      \def\svgwidth{.45\textwidth}
      \subfigure[]{\includegraphics[width=.35\textwidth,keepaspectratio]{figures/close/figclose1.pdf}}&
      \subfigure[]{\includegraphics[width=.35\textwidth,keepaspectratio]{figures/close/figclose2.pdf}}\\
    \end{tabular}
  \end{center}
  \caption{Case of \red{closely} situated objects. In this situation,
    it is impossible to fit object $i$ (in green) inside a rectangle
    $D_i$ that does not contain any points of objects $j \ne i$ (in
    red). A first idea would be to take a more complex contour for
    $D_i$ (e.g. the light green contour in $(a)$). A more convenient
    method is to add discretization points of the objects $j \ne i$ to
    $D_i$ , as shown in $(b)$.}
  \label{fig:extcontour}
\end{figure}

\begin{figure}
  \begin{center}
    \begin{tabular}{cc}
      \subfigure[]{\includegraphics[width=.45\textwidth,keepaspectratio]{figures/npts/fignptscontour.pdf}}&
      \subfigure[]{\includegraphics[width=.45\textwidth,keepaspectratio]{figures/npts/fignptsevolution.pdf}}\\
    \end{tabular}
  \end{center}
  \caption{Scattering problem with 64 objects. $(a)$: geometry, with initial discretization points in black and final skeleton points in red, $(b)$:
    evolution of the number of discretization points with the level in
    the hierarchical decomposition. The approximation error on a
    circle containing the scatterers is $10^{-7}$. Figures are generated by the {\sc Matlab} file
    \texttt{getnpts(64,1.)}.}
  \label{fig:64obj}
\end{figure}


\subsection{Size of the Approximate Scattering Matrix}

Let $N$ be the number of discretization points on each
contour, and let $N_0 = 2^n N $ be the total number of degrees of
freedoms (DOFs) in the initial discretization. At each level $0 \le l
\le n$ in the algorithm, there are $2^{n - l}$ objects with a total of
$N(l)$ DOF (i.e. $\approx 2^{-n + l}N(l) \red{equiv} \mathcal N(l)$
DOFs per object).

\red{The algorithm described above has been applied to different
  configurations.} For problems with $2^n$ closely located objects ($
0 \le n \le 8$) the total number of degrees of freedom $N(l)$ at each
iteration is displayed in figure~\ref{fig:effectn}. The initial number
of discretization points on each contour has been set to $N = 250$,
which is sufficient to resolve the length scales under considerations
here.  The \red{total} number of degrees of freedom is therefore $2^n
N$. When the frequency $k$ is of order 1 or lower (i.e. when the wave
length is of the order of the size of the object or larger) the final
number of points on the skeleton weakly depends on the initial number
of objects (see figures~\ref{fig:effectn}$(a,b)$), and it remains of
the order of $N$. As the frequency increases the final number of
skeleton points significantly increases: figure~\ref{fig:effectn}$(c)$
shows that more than $1000$ points are require\red{d} to represent the
scattering operator for $256$ objects with $k = 10$.

In algorithm~\ref{alg.multiple-scatterers}, the cost can be estimated
as follows:
\begin{itemize}
\item Assembly of the diagonal blocks $A_{ii}$ : $O(2^n N^2)$
  operations
\item At each level $l\ge 1$:
  \begin{itemize}
  \item building $2^{n - l+1}$ proxies
    \begin{itemize}
    \item Building $\v B$ and $\v C$: $O(2^{n-l} \mathcal N(l)M)$
    \item Building the interpolatory decomposition: 
    \item Computing $\v A^{-1}$ : $O(2^{n-l} \mathcal N(l)^3)$
    \item Building the proxy : $O(2^{n-l} \mathcal N(l)^2N(l+1))$
    \end{itemize}
  \item merging proxies into $2^{n - l}$ groups
    \begin{itemize}
    \item Building $\v A_{12}$ and $\v A_{21}$: $O(2^{n-l} \mathcal N(l)^2)$
    \item Merging blocks : $O(2^{n-l} \mathcal N(l)^3)$
    \end{itemize}
  \end{itemize}
\end{itemize}
so that to the leading order the total cost scales as
\begin{equation}
  \sum_l 2^{n-l} \mathcal N (l)^3 = \sum_l 2^{n-l} \lp
  \f{N(l)}{2^{n-l}}\rp^3  = \sum_l N(l)^3 2^{-2(n-l)}
\end{equation}
Figure~\ref{fig:effectn}$(d)$ represents this estimate as a function
of \red{$n$} for different frequencies. It shows that the present
method requires a number of operations that scales between $2^{3/2n}$
and $2^{2n}$.


\begin{figure}[t]
  \begin{center}
    \begin{tabular}{cc}
      \subfigure[$k = 0.1$]{\input{./figures/npts/effectn_k0.1_250.tex}}&
      \subfigure[$k = 1$]{\input{./figures/npts/effectn_k1_250.tex}}\\
      \subfigure[$k = 10$]{\input{./figures/npts/effectn_k10_250.tex}} &
      \subfigure[Estimated cost \red{for $2^n$ objects}]{\input{./figures/npts/costinv.tex}}\\
    \end{tabular}
  \end{center}
  \vspace{-0.3cm}
  \caption{$(a - c)$: Evolution of the number of degrees of freedom
    (DOFs) through the hierarchical solver for different frequencies
    $k$. The initial discretization (level 0) consists of $1$ to $256$
    objects with $250$ discretization points on each. Data generated
    with the {\sc Matlab} user-defined function
    \texttt{getnpts}. $(d)$: estimated cost taking into account only
    the computation of inverses, scaled by the estimated cost for one
    object. The solid black line corresponds to a cost proportional to
    $N_0$, the dashed black line corresponds to a cost proportional to
    $N_0^2$ and the dotted black line corresponds to a cost
    proportional to $N_0^{3/2}$. }
  \label{fig:effectn}
\end{figure}

\section{Concluding Remarks}
\label{sec.conclusion}

We have considered the scattering of a wave by solid smooth objects in
this work.  We show that the scattering matrix can be utilized to
accurately capture the relationship between the incoming and outgoing
wave outside a domain that encompass the scatterers. When the wave
length is large enough compared to the size of the scatterer, this
scattering matrix admits a low rank approximation, called the proxy
matrix, which allows it to represent the scattered field outside of
$D$ up to an arbitrary decomposition.  This proxy matrix can be
characterized in terms of the skeleton points, i.e. the set of points
of the initial discretization that are required to accurately
represent the far field. \red{The} influence of the different
parameters -- the number of discretization nodes, the frequency and
the shape of the object -- on the skeleton has been investigated, and
this approximation has been used in order to build a direct solver for
a large number of objects. A hierarchical method has been described in
which skeletons are build for each scatterer and then \red{merged} two
by two.

The method presented here is not fast, in the sense that the number of
operations required does not scale with the number of degrees of
freedom, but yet it allows the treatment of scattering by a large
number of objects with a reasonable cost. Furthermore, by not
requiring the full discretization matrix to be build at once, it
limits the memory requirements of the algorithm. Our ongoing efforts
are directed towards extending the present work to include scatterers
with corners.
\bibliographystyle{IEEEtran} % Include this if you use bibtex
\bibliography{bib/report-scatter-bib} % and a bib file to produce the
% end bibliography

\end{document}





